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A deterministic algorithm to compute 
approximate roots of polynomial systems 
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Pierre Lairez 


Abstract 

We describe a deterministic algorithm that computes an approximate root of n complex 
polynomial equations in n unknowns in average polynomial time with respect to the 
size of the input, in the Blum-Shub-Smale model with square root. It rests upon 
a derandomization of an algorithm of Beltran and Pardo and gives a deterministic 
affirmative answer to Smale’s 17th problem. The main idea is to make use of the 
randomness contained in the input itself. 


Introduction 

Shub and Smale provided an extensive theory of Newton’s iteration and homotopy continu¬ 
ation which aims at studying the complexity of computing approximate roots of complex 
polynomial systems of equations with as many unknowns as equations. 1 In their theory, an 
approximate root of a polynomial system refers to a point from which Newton’s iteration 
converges quadratic ally to an exact zero of the system—see Definition 1. This article answers 
with a deterministic algorithm the following question that they left open: 

Problem (Smale 2 ). Can a zero of n complex polynomial equations in n unknowns be found 
approximately, on the average, in polynomial time with a uniform algorithm? 

The term algorithm refers to a machine a la Blum-Shub-Smale 3 (BSS): a random access 
memory machine whose registers can store arbitrary real numbers, that can compute ele¬ 
mentary arithmetic operations in the real field at unit cost and that can branch according 
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to the sign of a given register. To avoid vain technical argumentation, I consider the BSS 
model extended with the possibility of computing the square root of a positive real number 
at unit cost. The wording uniform, algorithm emphasizes the requirement that a single finite 
machine should solve all the polynomial systems whatever the degree or the dimension. The 
complexity should be measured with respect to the size of the input, that is the number 
of real coefficients in a dense representation of the system to be solved. An important 
characteristic of a root of a polynomial system is its conditioning. Because of the feeling that 
approximating a root with arbitrarily large condition number requires arbitrarily many steps, 
the problem only asks for a complexity that is polynomial on the average when the input is 
supposed to be sampled from a certain probability distribution that we choose. The relevance 
of the average-case complexity is arguable, for the input distribution may not reflect actual 
inputs arising from applications. But yet, average-case complexity sets a mark with which 
any other result should be compared. 

The problem of solving polynomial systems is a matter of numerical analysis just as much 
as it is a matter of symbolic computation. Nevertheless, the reaches of these approaches 
differ in a fundamental way. In an exact setting, having one root of a generic polynomial 
system is having them all because of Galois’ indeterminacy, and it turns out that the number 
of solutions of a generic polynomial system is the product of the degrees of the equations, 
Bezout’s bound, and is not polynomially bounded by the number of coefficients in the input. 
This is why achieving a polynomial complexity is only possible in a numerical setting. 

The main numerical method to solve a polynomial system / is homotopy continuation. 
The principle is to start from another polynomial system g of which we know a root g 
and to move g toward f step by step while tracking all the way to / an approximate 
root of the deformed system by Newton’s iteration. The choice of the step size and the 
complexity of this procedure is well understood in terms of the condition number along the 
homotopy path. 4 Most of the theory so far is exposed in the book Condition. 5 The main 
difficulty is to choose the starting pair (g, if). Shub and Srnale 6 showed that there exists good 
starting pairs, and even many, for some measure, but without providing a way to compute 
them efficiently. Beltran and Pardo 7 discovered how to pick a starting pair at random and 
showed that, on average, this is a good choice. This led to a nondeterministic polynomial 
average-time algorithm which answers Smale’s question. Biirgisser and Cucker 8 performed a 
smoothed analysis of the Beltran-Pardo algorithm and described a deterministic algorithm 
with complexity N°^ oglogN \ where N is the input size. The question of the existence of a 
deterministic algorithm with polynomial average complexity it still considered open. 

This work provides, with Theorem 23, a complete deterministic answer to Smale’s problem, 
even though, as we will see, it enriches the theory of homotopy continuation itself only 
marginally. The answer is based on a derandomization of the nondeterministic Beltran and 

4 Beltran and Pardo, “Fast linear homotopy to find approximate zeros of polynomial systems”; Biirgisser and 
Cucker, “On a problem posed by Steve Smale”; Shub, “Complexity of Bezout’s theorem. VI. Geodesics in 
the condition (number) metric”. 

’Biirgisser and Cucker, Condition. 

6 Shub and Smale, “Complexity of Bezout’s theorem. V. Polynomial time”. 

7 Beltran and Pardo, “Fast linear homotopy to find approximate zeros of polynomial systems”, “Smale’s 17th 
problem: average polynomial time to compute affine and projective solutions”. 

8 Biirgisser and Cucker, “On a problem posed by Steve Smale”. 
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Pardo’s algorithm according to two basic observations. Firstly, an approximate root of a 
system / is also an approximate root of a slight perturbation of /. Therefore, to compute an 
approximate root of /, one can only consider the most significant digits of the coefficients of /. 
Secondly, the remaining least significant digits, or noise, of a continuous random variable are 
practically independent from the most significant digits and almost uniformly distributed. In 
the BSS model, where the input is given with infinite precision, this noise can be extracted 
and can be used in place of a genuine source of randomness. This answer shows that for 
Smale’s problem, the deterministic model and the nondeterministic are essentially equivalent: 
randomness is part of the question from its very formulation asking for an average analysis. 
It is worth noting that the idea that the input is subject to a random noise that does not 
affect the result is what makes the smoothed analysis of algorithms relevant. 9 Also, the study 
of the resolution of a system / given only the most significant digits of / is somewhat related 
to recent works in the setting of machines with finite precision. 10 

The derandomization proposed here is different in nature from the derandomization 
theorem BPPr = Pr, 11 which states that a decision problem that can be solved over the 
reals in polynomial time (worst-case complexity) with randomization and bounded error 
probability can also be solved deterministically in polynomial time. Contrary to this work, 
the derandomization theorem above relies on the ability of a BSS machine to hold arbitrary 
constants in its definition, even hardly computable ones or worse, not computable ones which 
may lead to unlikely statements. For example, one can decide the termination of Turing 
machines with a BSS machine insofar Chaitin’s f 1 constant is built in the machine. 

Acknowledgment I am very grateful to Peter Biirgisser for his help and constant support, 
and to Carlos Beltran for having carefully commented this work. I thank the two referees for 
their meticulous reading and their insightful suggestions. 
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1 The method of homotopy continuation 

This part exposes the principles of Newton’s iterations and homotopy continuation upon 
which rests Beltran and Pardo’s algorithm. It mostly contains known results and variations of 
known results that will be used in the next part ; notable novelties are the inequality relating 
the maximum of the condition number along a homotopy path by the integral of the cube of 
the condition number (Proposition 7) and a variant of Beltran and Pardo’s randomization 
procedure (Theorem 9). For Smale’s problem, the affine setting and the projective setting 
are known to be equivalent, 12 so we only focus on the latter. 


1.1 Approximate root 


Let n be a positive integer. The space <C n+1 is endowed with the usual Hermitian inner 
product. For d £ N, let Hd denote the vector space of homogeneous polynomials of degree d in 
the variables Xq, ..., x n . It is endowed with an Hermitian inner product, called Weyl’s inner 
product, for which the monomial basis is an orthogonal basis and ||a;g 0 • • • || 2 = ^ a yy - 
Let d i,..., d n be positive integers and let H denote Hd 1 x • • • x Hd n , the space of all systems of 
homogeneous equations in n+1 variables and of degree di,... ,d n . This space is endowed with 
the Hermitian inner product induced by the inner product of each factor. The dimension n 
and the d,’s are fixed throughout this article. Let D be the maximum of all di s and let N 
denote the complex dimension of ji, namely 


N = 



Elements of ji are polynomial systems to be solved, and 2 N is the input size. Note that 2 ^ N , 
n 2 ^ N and D ^ N . 

For every Hermitian space V, we endow the set §(V) of elements of norm 1 with the 
induced Riemannian metric d§: the distance between two points x,y £ §(V) is the angle 
between them, namely cosd§(x,y) = R e(x,y). The projective space P(V) is endowed with 
the quotient Riemannian metric dp defined by 

dp(M, [y]) d = min d§(x,Xy). 


An element of / £ Ti. is regarded as a homogeneous polynomial function C ra+1 —> C n . 
A root —or solution, or zero —of / is a point ( g P“ such that /(£) = 0. Let V be the 
solution variety {(/,£) £ % x P” | f(z) = 0}. For z £ <C ra+1 \ {0}, let d f(z) : C n+1 C n 
denote the differential of f at z. Let z 1 - be the orthogonal complement of C z in C n+1 . If 
the restriction df(z)\ z ± : z ± —> C n is invertible, we define the projective Newton operator Af, 
introduced by Shub, 13 by 

N(f,z) d = z - df(z)\~l(f(z)). 

It is clear that Af(f, A z) = XNf(f, z), so A f(f, —) defines a partial function P” —> P n . 


Definition 1. A point z £ P n is an approximate root of / if the sequence defined recursively 
by zo = z and Zk+i = A f(f, zyj) is well defined and if there exists (gP n such that /(£) = 0 

1 “Belt ran and Pardo, “Smale’s 17th problem: average polynomial time to compute affine and projective 
solutions”. 

13 Shub, “Some remarks on Bezout’s theorem and complexity theory”. 


(Symbols in the 
margin mark the 
place where they 
are defined.) 

n 


D 

N 


V 


A f{f,z) 


4 



and dp(zk,C) ^ 2 1 2 dp(z,f) for all k ^ 0. The point f is the associated root of z and we 
say that z approximates f as a root of f. 

For / G PL and z G C " +1 \ {0}, we consider the linear map 

E(/,z) : (ui ,...,u n ) G C n >-+df(z)\~l ^\/di||^|| dl_ 1 ui,• • •, \fdn\\z\\ dn ~ 1 u r ^j G z 1 - 

and the condition number 1 ^ of f at z is defined to be p(f, z) d = ||/|| ||S(/,z)||, where \\^{f,z)\\ p(f,z) 
is the operator norm. When df(z)\ z ± is not invertible, we set p(f,z) = oo. The condition 
number is often denoted /i norrrl but we stick here to the shorter notation p. For all u, v G C x 
we check that p(uf,vz) = p(f,z). We note also that p(f,z ) ^ y/n ^ l . 15 The projective 
/x-theorein (a weaker form of the better known projective 7 -theorem) relates the condition 
number and the notion of approximate root: 

Theorem 2 (Shub, Srnale 16 ). For any (/, C) € V and z G P", if D 3 / 2 p(f, f)dp(z, f) ^ 
then z is an approximate root of f with associated root f. 

Remark. The classical form of the result , 17 requires D 3 / 2 p(f, f) tan(dp( 2 , f)) < 3 — y/7. 

The hypothesis required here is stronger: since D 3 / 2 p(f,f) ^ 1, if D 3 / 2 p(f,f)dp(z,f) ^ 5 
then dp(z,f) ^ 5 and then tan(dp(z, f)) ^ 3tan(|)dp(,j, f) < oW/pj c) because tan(f) 

3 — \fl. The symbol indicates an inequality that is easily checked using a calculator. 

The algorithmic use of the condition number heavily relies on this explicit Lipschitz 
estimate: 

Proposition 3 (Shub 18 ). Let 0 ^ e ^ f. For any f,g G P (PL) and x,y G P", if 
/j(/,s)max (p 1,2 d v (f 1 g) 1 D 3/2 d v (x,y) S j ^ | 
then (1 + e)~ V(/, x) < p(g, y) < (1 + e)p(f, x). 

1.2 Homotopy continuation algorithm 

Let JcKbe an interval containing 0 and let f G / H> f t G P (PL) be a continuous function. 

Let f be a root of fo such that d/o(C)|c J - is invertible. There is a subinterval J C I containing 0 
and open in I, and a continuous function f G J 4 (t G P" such that Co = C and ft(ft) = 0 
for all t G J. We choose J to be the largest such interval. 

Lemma 4. If t ^ ft is C 1 on I and if p(ft, ft) is bounded on J, then J = I. 

Proof. Without loss of generality, we may assume that I is compact, so that ||/ t || is bounded 
on I. Let M be the supremum of p(ft, Ct)ll/t|| 011 J- From the construction of ft with the 
implicit function theorem we see that i G f 4 f t is M-Lipschitz continuous. Hence the 
map t G J H > ft extends to a continuous map on J. Thus J is closed in J, and I = J 
because J is also open. □ 

14 Shub and Smale, “Complexity of Bezout’s theorem. I. Geometric aspects”; see also Biirgisser and Cucker, 

Condition , §16, for more details about the condition number. 
lo Biirgisser and Cucker, “On a problem posed by Steve Smale”, Lemma 16.44. 

16 Shub and Smale, “Complexity of Bezout’s theorem. I. Geometric aspects”. 

17 Blum, Cucker, Shub, and Smale, Complexity and real computation, §14, Theorems 1 and 2. 

18 Shub, “Complexity of Bezout’s theorem. VI. Geodesics in the condition (number) metric”, Theorem 1; see 
also Biirgisser and Cucker, Condition, Theorem 16.2. 
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Proposition 5. Let (/, C) P {%) and 0 < e < If D 3 / 2 p(f, £) 2 d P (/, g) s£ 

then: 

(i) there exists a unique root g of g such that dp(£, 77 ) ^ (1 + e)p(f, £)dp(/, g); 

(a) (1 + e)-V(/, C) < v) < (1 + e)M/> 0 ; 

(Hi) £ approximates g as a root of g and 77 approximates ( as a root of f. 

Proof. Let t £ [0, 1] f t € P (H) be a geodesic path such that /o = f, fi = g and ||/ t || = 

dp(f, g)- Let t £ J 1 —> (t be the homotopy continuation associated to this path starting from 
the root £ and defined as above on a maximal interval J C [0,1], Let pt denote p{ft, (t)- 
For all t £ J we know that ||£t|| ht\\ft\\, 19 so that 

\\du ^ d r (f,g) / p u du. ( 1 ) 

Jo 

Let J' be the closed subinterval of J defined by J' = [t £ J | Mt' < t, D 3 ! 2 /i 0 dp(£o, (t>) < |}. 

For all t £ J' we have D 3 / 2 /r 0 dp(Co, £t) < f, by definition, and D 1 / 2 p 0 d P (f 0 , f t ) < D 3 / 2 p( j d P (f, g) < 

|, by hypothesis. Thus, Proposition 3 ensures that 

(1 + e)“ Vo < Pt < (1 + e)A*o, for all t £ J'. (2) 

Thanks to Inequality (1) we conclude that dp(£o,£t) ^ (1 + e)t dp(f, g)po, for all t £ J\ so 
that D 3/ V 0 dp(Co,Ct) < 7 , usin S the assumption D 3 / 2 p(fX) 2 d P (f,g) ^ This proves 

that J' is open in J. Since it is also closed, we have J' = J. Since p t is bounded on J', by 

Inequality (2), Lemma 4 implies that J' = J = [0,1]. Now, Inequalities (1) and (2) imply 
that dp(£o,£i) ^ (1 + e)d P (f, g)po- This proves (i) and (ii) follows from (2) for t = 1. 

To prove that 77 approximates £ as a root of /, it is enough to check that 

D 3 / 2 p{f,Qd P {C,v) < (l + e)-D 3/ V(/V)%(/,sO < | ^ 

by Theorem 2. To prove that £ approximates 77 as a root of < 7 , we check that 

D 3/2 p(g, ??)dp(£, 77 ) s£ (l + £) 2 D 3/ V(/,£) 2 d P (/,V < pt- 

This proves (iii) and the lemma. □ 

Throughout this article, let£=pp^. = V>-S = Vy and B' = -J r. The main result that A, B, B ', £ 
allows computing a homotopy continuation with discrete jumps is the following: 

Lemma 6 . For any (/, £) £ V and g £ H and for any z £ P™, if D 3 ! 2 p(f, z)d P (z, £) < A 
and D 3 / 2 p(f, z) 2 d P (f, g) < B' then: 

(i) z is an approximate root of g with some associated root g; 

(ii) (1 + e)~ 2 p(f , z) < p(g, 77) < (1 + e) 2 p(f, z); 

(iii) D 3 / 2 p(g,g)d P (z,g) 

If moreover D 3 / 2 p(f, z) 2 d P (f,g) < B iften: 

19 Biirgisser and Cucker, Condition, Corollary 16.14 and Inequality (16.12). 
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(iv) D 3 / 2 g{g, z')dp(z', g) ^ A, where z' = J\f(g,z). 


Proof. Firstly, we bound g(f,(). Since D 3 / 2 g(f, z)dp(z, () < A = |, Proposition 3 gives 

(1 + e) _ V(/, C) < z ) < (! + e)g(f, C). 


Next, we have D 3 / 2 g(f, Q 2 dp(f, g) ^ (1 + £) 2 iT , thus Proposition 5 applies 

and £ is an approximate root of g with some associated root g such that dp(f,g) ^ (1 + 
e)h(fX)dp(f,g) and (1 + e)~V(/, C) < h{g,v) < (1 +e)/x(/,C) and this g ives (“)• 

Then, we check that 2 approximates as a root of g. Indeed 


dr{z, g) < d r (z, Q + d r (C, rj) < 


A + (1 + s) 2 B' 
D 3 / 2 g(f,z) 




(1+£) 2 {A + {1 + £) 2 B') 

D 3/2 g{g,g) 


And (1 + e) 2 (A + (1 + £) 2 B') A < |, so Theorem 2 applies and we obtain (i) and (iii). 

We assume now that D 3 / 2 g{f, z) 2 dp(f,g) ^ B. All the inequalities above are valid with B’ 
replaced by B. By definition of an approximate root dp(z',g) ^ \dp{z,g), so that 

D 3/2 g{g,g)d P (z',g) < ^(1 + e) 2 (A + (1 + £) 2 B) 

Thus (1 + e)~ l g{g, rj) < g(g, z') < (1 + e)g{g, g). 

To conclude, we have D 3 / 2 g(g,z')d(z , ,g) < |(1 +e) 3 (A + (1 + £) 2 B) A. □ 


Let f,g € S (ft), with / 7^ —g. Let t G [0,1] >->• T(g, /, t) be the geodesic path from g to / 
in S(fL). The condition f —g guarantees that the geodesic path is uniquely determined, 
namely 


r ( g ./,t) = 3m((1 :' t)a) g + ^/. 


(3) 


sin(a) sin(a) 

where a = d§(f,g) G [0, n) is the angle between / and g. 

Let z G P” such that D 3 ^ 2 g(g, z)dp(z,g) ^ A, for some root g of g. By Lemma 6(i), 
applied with g = / and g = £, the point z is an approximate root of g , with associated root g. 
Given g and z, we can compute an approximate root of / in the following way. Let go = g, 
to = 0 and by induction on k we define 


hk — hifk, Z k ), tk-(-1 — tk T 


B 

D 3 / 2 g,^d§(f, g) ’ 


gk+ 1 = T(g,f,t k+1 ) and z k + i = A/^fc+i, Zfc)- 


Let K(f, g, z), or simply K , be the least integer such that tj-c+i > 1, if any, and K(f, g, z) = oo 
otherwise. Let M(f,g,z) denote the maximum of all g k with 0 ^ k ^ K. Let HC be the 
procedure that takes as input /, g and z and outputs zk- Algorithm 1 recapitulates the 
definition. It terminates if and only if K < oo, in which case K is the number of iterations. For 
simplicity, we assume that we can compute exactly the square root function, the trigonometric 
functions and the operator norm required for the computation of g(f, z). Section 2.5 shows 
how to implement things in the BSS model extended with the square root only. 

Let ht = T(f,g,t) and let t G J K»• ft be the homotopy continuation associated to f G 
[0,1] i— y ht, where go is the associated root of z, defined on a maximal subinterval J C [0,1], 
Let 

M{f,g,z) = f sup g(ft, ( t ) and I p (f,g,z) d = [ g(h t , g t ) p dt. (4) 

teJ J j 


r (g,f,t) 


K{f,g,z) 
M(f,g,z) 
HC (f,g,z) 


M(f,g,z) 


7 



Algorithm 1. Homotopy continuation 
Input. /, g £ SfH) and 2 £ P™. 

Precondition. There exists a root g of g such that 52 D 3 l 2 p(g, z)dp(z, rj) ^ 1. 

Output, w £ P n 

Postcondition, w is an approximate root of /. 
function HC(/, g, z) 

1 / (l0in 3 / 2 //(g,z) 2 d s (/,g)) 

while 1 > t do 

h <- r (g,f,t) 
z «— A/"(/i, z) 

t <- t + 1/ (l01L> 3/2 /r(/i, z) 2 d$(f , g)) 

end while 
return z 
end function 


The behavior of the procedure HC can be controlled in terms of the integrals I p (f, g, z). It 
is one of the corner stone of the complexity theory of homotopy continuation methods. The 
following estimation of the maximum of the condition number, along a homotopy path, in 
terms of the third moment of the condition number seems to be original. It will be important 
for the average complexity analysis. 


Proposition 7. If J = [0,1] then M(f,g,z) ^ 151 D 3 / 2 I 3 (f,g, z). 


Proof. Let e = j and let s £ [0,1] such that g(f s ,(s) is maximal. For all t £ [0,1], 
<k{fs,ft) < 1 1- s\d§(f,g). Thus, if 

^ 4(l + £)D 3 / 2 /i(/ s ,C s ) 2 ds(/,ff)’ ^ 

then n(ft, Ct) ^ (1 + e)~ 1 p(f s , CO, by Proposition 5. Since d§(f,g) < 7 t, the diameter of the 
interval H of all t £ [0,1] satisfying Inequality (5) is at least 47 r ( 1 +g )j)3 e /V(/ ^ y Thus 


[ M/t,Ct) 3 d£ ^ / 

Jo Jh 


v{fs,Cs)' c 

(l + e ) 3 


d t ^ 


ep(fs, C) 




1 M/»> C) 


4tt(1 + £ ) 4 L» 3 /2 151 D 3 / 2 


□ 


Theorem 8 (Shub 20 ). With the notations above, if D 3 ^ 2 p(g,z)dp(z,r]) ^ A then: 


(i) HC (f,g,z) terminates if and only if ^(f, g, z) is finite, in which case J = [0,1]; 
If moreover HC(/, g, z) terminates then: 

(ii) (1 + £)- 2 M(f,g,z) < M(f,g,z) ^ (1 + e) 2 M(f,g,z). 

(in) K(f,g,z) ^ 136 D 3 / 2 ds{f,g)I 2 (f,g,z); 


(iv) HC(/, g, z) is an approximate root of f; 

(v) D 3 / 2 p.{fX)d P (RC(f,g,z),C) ^ where ( is the associated root o/HC(/, g, z). 


20 Shub, “Complexity of Bezout’s theorem. VI. Geodesics in the condition (number) metric”. 



Proof. Let T] k denote f tk - Since D 3 / 2 n k dp(g k , gk+i) =SS B for all k ^ 0, Lemma 6 (iv) proves, 
by induction on k that D 3 ^ 2 g k dp(z k , g k ) ^ A for any k ^ 0. 

Assume that [0, tfc] C J for some k ^ 0 and let t £ [t k ,t k +i] H J so that 

D 3/2 g 2 k d(g k ,h t ) ^ D 3/2 g 2 k d(g k ,g k+1 ) ^ B. 

Because D 3 / 2 g k d(z k , g k ) ^ A, Lemma 6 (ii) applies to (g k . r/ k ), ht and z k and asserts that 

(1 + e)~ 2 g k ^ g(h t , C t) < (1 + e) 2 g k . (6) 


By definition g k (t k + 1 — t k ) = fl3/2 ^^ , so integrating over t leads to 
[ n{h t ,CtfAt ^ (l + ff)“ 4 ^/r|(^ + i -tj) = 

j 0 ;_n 


3=0 


- h) - 2 ±$ 


and 


(1 + e) 4 D 3 / 2 d§(/, 5 )’ 

(l + e) 4 (fc + l )B 


if, 9) 


(7) 

( 8 ) 


Assume now that / 2 (/, g, z) is finite. The left-hand side of Inequality (7) is finite so there 
exists a k such that t k +i fL J. But then Inequalities ( 6 ) shows that /i ( is bounded on J which 
implies, Lemma 4 that J = [0,1]. And since t k +i ^ J , this proves that K is finite. 

Conversely, assume that K is finite, i.e. HC (f,g,z) terminates. Then there exists a 
maximal k such that [0, t k \ C J and thus for all t £ J 


d(htXt) < (1 + s) 2 ma xn(g k , z k ). 


So g(ht,Ct) is bounded on J, which implies that J = [0,1], and thus k = K. Inequality ( 8 ) 
then shows that I 2 (f,g,z) is finite, which concludes the proof of (i). We keep assuming 
that I\ is finite. Inequality ( 6 ) shows (ii). Since [0,t^] C [0,1], by definition, Inequalities (7) 
and ( 8 ) shows that 

B(l \ £ y D 3 / 2 ds(f,g)h(f,g,z) - 1 < K < ^^D 3 / 2 d s (f,g)I 2 (f,g,z). 

We check that ^ 136, which gives (iii). Finally, Lemmas 6 (i) and 6 (iii) show that zk 

approximates Cl as a root of / and that D 3 t 2 g{f, (i)dp(zK, Ci) ^ 2 3 ; which gives (iv) 
and (v). □ 


1.3 A variant of Beltran-Pardo randomization 

An important discovery of Beltran and Pardo is a procedure to pick a random system and 
one of its root simultaneously without actually solving any polynomial system. And from 
the complexity point of view, it turns out that a random pair ( g , 77 ) £ V is a good starting 
point to perform the homotopy continuation. 

Let g £ S("H) be a uniform random variable, where the uniform measure is relative to the 
Riemannian metric on S(77). Almost surely g has finitely many roots in P". Let 77 be one of 
them, randomly chosen with the uniform distribution. The probability distribution of the 
random variable (< 7 , 77 ) £ V is denoted p s td- The purpose of Beltran and Pardo’s procedure 21 

- 1 Be]Iran and Pardo, “Fast linear homotopy to find approximate zeros of polynomial systems”, §2.3; see also 


Pstd 
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is to generate random pairs (g,rf), according to the distribution p st d without solving any 
polynomial system. We give here a variant which requires only a uniform random variable 
in S(C W ) ~ S("H) as the source of randomness. 

Let us assume that / = (/i,..., /„) £ S("H) is a uniform random variable and write / as 

n 

fi — T fdi.X q ^ ' Cli,jXi -p f,j (xo, • . • , X n ), 

3 = 1 


for some c^, oqj £ C and /' G H r [ t such that f-(e o) = 0 and d/'(eo) 
(/(,..., f n ) £ H. By construction, f(e o) = 0 and d/'(eo) = 0. Let 


M d = 


( ai 7 i 

®l,n 

\^n,l 

^n,n 


Cl 


(p>nx (n+1) 


0. Let /' = 


Almost surely, ker M has dimension 1; Let ( £ P" be the point representing ker M and 
let f £ S(C ra+1 ) be the unique element of ker M nS(C" +1 ) whose first nonzero coordinate is 
a real positive number. Let = ('Ll, • • •,T,,) £ T~L be defined by 


'I',: 


di -1 


‘= \[di ( '}Zx t C t 


777-0 n X q 


\i —0 


3=0 


(9) 


where Q denotes the complex conjugation. By construction ’J / m,c , (C) = 0- Let u £ U(n+ 1 ), 
the unitary group of C n+1 , such that u{e o) = (. The choice of u is arbitrary but should 
depend only on £. For example, we can choose u, for almost all £, to be the unique element 
of U(n + 1) with determinant 1 that is the identity on the orthogonal complement of {eo, C} 
and that sends eo to (. Finally, let g = f o u -1 + T m c G 'H- By construction g(Q = 0. We 

def ,S 

define BP(/) = (g, which is a point in the solution variety V. 


Theorem 9. If f £ S (H) is a uniform random variable, then BP(/) ~ p s td- 


Proof. We reduce to another variant of Beltran-Pardo procedure given by Biirgisser and 
Cucker 22 in the case of Gaussian distributions. Let / € §(7f) be a uniform random variable, 
and let \ £ [0, oo) be an independent random variable following the chi distribution with 2N 
degrees of freedom, so that xf is a centered Gaussian variable in H with covariance matrix I 2 N 
(which we call hereafter a standard normal variable). For £ £ P”, let C P be the subspace 
of all g such that g(Q = 0 and dg(£) = 0 and let be the orthogonal complement of Rq 
in H. The system xf splits orthogonally as xf + X^, where xf G R e0 an d € S eo are 
independent standard normal variables. 

Let M £ C nx ( n+1 \ ( £ P", f £ § n and u £ U(n+ 1) be defined in the same way as in the 
definition of BP(/). The map that gives M as a function of h is an isometry S eo —> C rax ( n + 1 ) ) 
so x'Af is a standard normal variable that is independent from /'. Let A £ S(C) be an 
independent uniform random variable, so that Af is uniformly distributed in ker M (~l S™ 
when M has full rank, which is the case with probability 1. The composition map g £ R eo <—> 
g o u £ R^ is an isometry. Thus, conditionally to (, the system xf 0 W 1 is a standard 

Biirgisser and Cucker, Condition , Chap. 17. 

22 Biirgisser and Cucker, Condition , Theorem 17.21(a). 


BP(/) 
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normal variable in R q. As a consequence, and according to Brirgisser and Cucker 23 , the 
system 

H = f xf ° w -1 + ^ x m,\C 

is a standard normal variable in R and £ is uniformly distributed among its roots. Hence i?/||f?| 
is uniformly distributed in § (R) and (H/\\H\\, f) ~ p st d- 

We check easily that ||1E , m,ac'II = ||M||i? = ||/i||, where ||M||i? denotes the Froebenius 
matrix norm, that is the usual Hermitian norm on C IX (” +1 ). Moreover \\f o u -1 || = ||/'||, 
this is the fundamental property of Weyl’s inner product on H. Thus ||fJ|| = ||/|| = y, and 
in turn 

(/' ° u- 1 + 'kM.AC', C) = (H/\\H\\, C) ~ Pstd, 

which is almost what we want, up to the presence of A. Let A £ C nxn be the diagonal matrix 
given by It is clear that Tm.ac' = The map M H>• AM is an isometry 

of <C rax ( n+1 ) and ker AM = kerM so and (yA M,u,Cf) have the same probability 

distribution. Since x..f is independent from yM and A, it follows that the system H ' defined 

by 

has the same probability distribution as H. To conclude the proof, we note that ||if'|| = y 

and that {H'/x, C) = BP(/). □ 


Given / £ S('H), Beltran and Pardo’s algorithm proceeds in sampling a system g £ §(H) 
from the uniform distribution and then computing HC(/, BP(g)). If the input / is a uniform 
random variable then we can evaluate the expected number of homotopy steps E(if(/, BP(g))). 
Indeed, let g be root of g, uniformly chosen, the theorem above asserts that BP(g) has 
the same probability distribution as (g,rf) so E(A"(/, BP(<?))) = E (K(f,g,g)). Thanks to 
Theorem 8(iii), it is not difficult to see that E (K(f,g,rj)) ^ 214 D 3 t 2 E(p(g, g) 2 ). This is why 
the estimation of E(p,(g,g) 2 ) is another corner stone of the average complexity analysis of 
homotopy methods. Deriving from a identity of Beltran and Pardo, we obtain the following: 


Theorem 10. If (g,g) ~ p 8 td, then E(p(g,g) p ) ^ nN) p/2 for any 2 < p < 4. 

Proof. Let s = p /2 — 1. Beltran and Pardo 24 state that 




r(iv - s) 


k —1 


+ 1\ Tfk — s) k+s 


k + 1 


) m 


We use the inequalities x y r(;r) ^ T(a; — y) < (x — 1) y r(:r), for x £ [1, oo) and y £ [0,1], 
which comes from the log-convexity of T. In particular, since 0 ^ s < 1, 


r(iv-s) r(iv-s) 


Thus 


E(y(g,g) 2+2s ) < N 1+s 


■ + 1 W -s).-i 

r(i) 


m 


~ s n~ k+s 


23 Biirgisser and Cucker, Condition , Theorem 17.21(a). 

24 Beltran and Pardo, “Fast linear homotopy to find approximate zeros of polynomial systems”, Theorem 23. 
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On the one hand (1 — s)T(l — s) = F(2 — s) ^ T(2) = F(l), so 


n+ h r ( i - , ) n - 1 « ( " +1)n -2_ n — < - 


1+s 


2 J r(i) 


2 1 - s 


On the other hand, 


j:( n k + + y k - ir ' n ^ + ' in ' + ’ s( n r' 1 ”" 

k—2 ' ' ' k —3 ' 


n+1 


= n 1+s 1 + 


1 \ 


n+1 


1 — 5 


- 1 - 


n + 1 1 /n + P 

n n 2 V 2 , 


-:— + 


4 4(1 — s) 

Putting together all above, we obtain the claim. 


□ 


2 Derandomization of the Beltran-Pardo algorithm 

2.1 Duplication of the uniform distribution on the sphere 

An important argument of the construction is the ability to produce approximations of two 
independent uniform random variables in S 2W_1 from a single uniform random variable 
in S 2Ar - x given with infinite precision. More precisely, let Q be a positive integer. This section 
is dedicated to the construction of two functions |_—J q and {—}q from the sphere S 2W_1 to 
itself, respectively called the truncation and the fractional part at precision Q. For u £ S 2Ar_1 , 
[uJq is close to u and if u is uniform random variable, then {u}q is nearly uniformly 
distributed in § 2Ar_1 and nearly independent from |+Jq> hi the following sense: 

Lemma 11. For any u £ S 2Ar_1 , d§ (LuJq,u) ^ 37V 1 / 2 /Q. Moreover, for any continuous 
nonnegative function 0 : S 2Ar_1 x S 2Ar_1 — > R, 


1 


voi(s 2Ar_i ) y§ 2 jv- 


1 0 (Wq.Wq) du 


< 


exp(^ 
vol(S 2iV 


-n 2 


1§2N-1 J§2N-1 

def 


0 (|_uj q , u) dudu. 


For i£l, let A(x) denote the integral part of a and let Aq(o) =‘ Q 1 A(Qa) be the trunca¬ 
tion at precision Q. For x £ R 2Ar + let Aq(x ) £ R 2 ^ 1 be the vector (Aq(x i), ..., Aq(x2n-i)) 


and let Bq(x) = f (x — Aq(x))Q, which is a vector in [0, 1] 2JV_1 . We note that \\Aq{x) — x \\ 2 < 
(27V — 1)/Q 2 , because the difference is bounded componentwise by 1 /Q. 

Let C and C + denote [—1, l) 2Ar_1 and [0, l) 2Ar_1 respectively, and let F{x) = (l + ||:r|| 2 ) _Ar . 
We first show that if x £ C is a random variable with probability density function F (divided 
by the appropriate normalization constant) then Bq(x) is nearly uniformly distributed in C+ 
and nearly independent from Aq(x). 

Lemma 12. For any continuous nonnegative function 0 : [— 1, l]™" 1 x [0, l]™" 1 -+ R, 

J e(A Q (x),B Q (x))F(x) dz^exp^^^) J J Q(A Q (x),y)F(x)dxdy. 
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Proof. For any integers — Q ^ ki < Q, for 1 ^ ^ 27V — 1, the function Aq is constant 


X, 


(2 Q ) 2 


on the set Y\a=i 1 Qd > an d these sets form a partition of C. Let Xi,..., 
denote an enumeration of these sets and let a k denote the unique value of Aq on X k . The 
diameter of Xk is y/2N — 1 /Q. Since the function x G [0, oo) i-A —N log(l+a: 2 ) is 7V-Lipschitz 
continuous, we derive that 


maxF < e NV 2 N~i/Q minF e 2 JV 3/2 /Q m [ n p. 

Xk Xk Xk 


( 10 ) 


For any 1 ^ k ^ (2 Q) 2N 1 , we have 

0 (Aq(x),B q (x)) F(x)dx ^ maxF / 0 (a k , {x - a k )Q) dx, 

X k JXu 


>x k 


because Aq{x ) = a k on Xk and by definition of Bq(x). A simple change of variable shows 
that 

/ Q(a k ,(x-a k )Q) dx = vol(X k ) Q(a k ,y)dy. 

Jx k Jc + 

Besides, for all y G C + , 

1 


Q(a kl y) < 


0 {A Q {x),y) F{x)dx. 


vol(A fe ) min A y F J Xk 
Putting together all above and summing over k gives the claim. 


□ 


Thanks to a method due to Sibuya, we may transform a uniform random variable of (7+ 
into a uniform random variable in § 2JV_1 . Let x — (aq,..., aqjv-i) G < 7 - 1 -, let yi,... , y X - 1 
denote the numbers xjv+i> ..., X 2 N -1 arranged in ascending order, and let = 0 and y x = 1 . 

Let S{x) G R 2W denote the vector such that for any 1 ^ i ^ N S(x) 

S(x) 2 i -i = f \Jyi ~ Vi-i cos( 2 TTXi) and S(x) 2i = f \Jyi ~ Vi -1 sin(27nEj). ( 11 ) 

Proposition 13 (Sibuya 25 ). If x is a uniformly distributed random, variable in C + , then S(x) 
is uniformly distributed in S 2W_1 . 

We now define |_—J q and {—}q• Let LI C R 2JV be the set of all x G R 2Af such that Ha^oo = 1. 

It is divided into 47V faces that are isometric to C: they are the sets Ef = {x G E | x* = e}, 
for £ G { — 1,1} and 1 ^ i ^ 27V and the isometries are given by the maps 

• E i y (7, X I ^ (x\, . . . , Xi— i, Xi- (_i, - - - , X n ). 

Through these isometries, we define the functions A'q and B'q on E: for x G Ef we set 

Aq(x) = t- l £ {A Q {t^{x))) G Ef and B' Q {x) d = B Q (U, e (x)) G C+. 

Let isqo : u G S 2A,_1 i-A u/HuHoo G S and its inverse v 2 : x G E a?/||a;|| G § 2JV_1 . Finally, 

we define, for u G S 2 ^' 1 , using Sibuya’s function S, see Equation (11), L u Jq> { u }q 

\u\ Q d = v 2 (A' q {v 00 {u))) and {u} Q d = S (Bq^^u))) . ( 12 ) 

We may now prove Lemma 11. 


25 Sibuya, “A method for generating uniformly distributed points on A-dimensional spheres”. 
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Proof of Lemma 11. Let u £ S 2JV 1 . It is well-known that d§(|_itjQ, w) < §||[wJq — u \\- 
Furthermore, the map i / 2 is clearly 1-Lipschitz continuous on E so || |_wJq — u|| ^ ll^g( iy oo(w)) — 
foo(u)\\ and we already remarked that the latter is at most \/2N — 1 /Q. With | \[2 3, 

this gives the first claim. 

Concerning the second claim, we consider the partition of the sphere by the sets iG>(Ef). 
For any u £ ^ 2 (EIf), we have 


^oo(w) 


( Mi— 1 'U'i+l Un 

I , . . . , , 1, , . . . , 

V Ui Ui Ui Ui 


Hence, by Lemma 14 below, the absolute value of the Jacobian of the map : zq(E?) —> Ef 
at some z/ 2 ( x ), x € Ef, is precisely 

{y 2 {x)i)~ 2N = ||a:|| 2Ar = (1 + ||t i) e(a;)|| 2 ) Ar = F(t ite (x))~ l . 


Thus, the change of variable gives 



0 i v 2 { a 'q{x)) , S(B' q (x))) F(ti, e (x))dx, 


and then Lemma 12 (applied over Ef with the isometry : Ef —> C) implies 


^ exp 



0 (f 2 ( a q( x )) > S(y)) dy F(t i>e (x))dx 


and Proposition 13 gives the equality 



vo^S 2 ^ 1 ) 


/ s 2iV - 


0 { v 2 ( a q( x )) > v ) dv F(t i>e (x))dx, 


and applying the inverse change of variable '■ Ef —> i/ 2 (Ef) and summing over i and e gives 
the claim. □ 


Lemma 14. Let H = {a: £ R 2N | x\ > 0} and let p> be the map 

p : S 2Ar_1 n H —» K 2iv_1 



Then, for any u £ S 2Ar 1 fl H , | det (dc^(xt)) | = iq 2N , where E> 2N 1 and R 2N 1 are considered 
with their usual Riemannian structures. 


Proof. Let if : u £ R 2N fl H — > ( qyyj, so that tp is the restriction of 'if to the 

sphere § 2Ar_1 . Firstly, the matrix of d if{x), for some x £ R 2Af , in the standard bases of M 2JV 
and R 2N ~\ is given by 


Mat (d</?(a;)) 


/ ~X2 


Xi 


\-X2N o 



(13) 


Let u G S 2Ar 1 D H. We may assume without loss of generality that u is of the form (iq, w 2 , 0,..., 0), 
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with u\ + U 2 = 1, because |det(dy>(u))| is invariant under any unitary transformation of u 
that preserves the first coordinate. Let T C K 2JV be the tangent space at u of S 2Ar ~b 
Naturally, d<p(rt) = df/’(u)|T- An orthonormal basis of T is given by {/, e ^,..., e 2 jv}, 
where / = (— 112 , u 1 ,0,... ,0) and where ei is the ith coordinate vector. Using Equation (13), 
we compute that d= uf 2 e 1 and that d <p(u)(ei) = for 3 ^ i ^ 2N. 

Thus |det (di^(u))| = uf 2N . □ 

The orthogonal monomial basis of H gives an identification PL ~ R 2Ar and we define this 
way the truncation |_/J q and the fractional part {f} q of a polynomial system / € S (PL). 
The derandomization relies on finding a approximate root of |_/Jq, for some Q large enough, 
and using {/}q as the source of randomness for the Beltran-Pardo procedure. Namely, we 
compute HC([/Jq, BP({/}g)). Almost surely, this computation produces an approximate 
root of L/Jq- If Q is large enough, it is also an approximate root of /. The main technical 
difficulty is to choose a precision and to ensure that the result is correct while keeping the 
complexity under control. 

2.2 Homotopy continuation with precision check 

Let /, /', g £ S(PL) and let 77 £ P n be a root of g. Throughout this section, we assume that 
ds(/, f) ^ p, for some p > 0 and that d§(f, g) ^ 7 r/ 2 . Up to changing g into — g , the latter is 
always true, since d§(f, —g) = n — d§(f,g). The notations I 2 , M and M used in this section 
have been introduced in §1.2. If p is small enough, then Y{C{f',g 1 ri) is an approximate root 
not only of /' but also of /. But if p fails to be small enough, HC(/ , ,< 7 ,??) may not even 
terminate or, to say the least, HC(/ , ,g, 77 ) may take arbitrarily long to compute something 
that is not an approximate root of /. To control the complexity of the new algorithm, 
it is important to be able to recognize this situation at least as fast as HC(/, g, rf) would 
terminate. 

As in §1.2, let f t = T{g 1 f 1 t) and f[ — r(g,f,t). Let t £ J —> Q. £ P n be the homotopy 
continuation associated to f t , on [0,1], and t £ J' —> (,[ £ P" be the one associated 
to fj-, defined on some maximal intervals J, J 1 C [0,1] containing 0. Let p, t = 
and p! t = Q. 

Lemma 15. d s (f t , /() < 2 d§(f, f) for any t £ [0,1]. 

Proof. Let a t = d§(f t , f[), f3 = d§(f, g) ^ | and 7 = dg(/', g). Without loss of generality, we 
may assume that ay < f, otherwise the inequality a t ^ 2 ai that we want to check is trivial. 
The spherical law of cosines applied to the spherical triangle {g, f t , //} gives the equality 

cos a t = cos(f/3) cos(t 7 ) + sin(t/?) sin(t 7 ) cos A, (14) 

where A is the angle of the triangle at g. We deal with three cases. Firstly, we assume 
that 7 ^ (p Then cosa t decreases at t increases: Indeed, Equation (14) rewrites as 

cosa t = cos (t/3 — t'y) — sin(t/3) sin(f 7 )(l — cos A) (15) 

and, as t increases, cos (t/3 — t'y) decreases, because |/3 — 7 | ^ 7 r, and both sin(t/3) and sin(f 7 ) 
increase, because /3 ,7 ^ |. Thus cosa t ^ cosai, for 0 ^ t < 1, and it follows that a t < ai- 
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Second case, we assume that 7 > f and /3 = f. For t £ [0,1], Equation (15) shows that 

cosa t ^ cos(| — 7 ) — (1 — cos A) = sin 7 + cos A — 1, 

using cos(f/3 — 1 7 ) ^ cos(/3 — 7 ) and 1 — cos(A) ^ 0. Equation (14) shows that cosoq = 
sin 7 cos A. In particular cos A ^ 0, since oq ^ f and sin 7 ^ 0. It follows that 

2 sin 2 7 cos 2 A ^ sin 4 7 + cos 4 A ^ sin 7 + cos A, 

and finally that cos(2oq) ^ cosaq, because cos(2aq) = 2cos 2 oq — 1. Since 2aii ^ 7 r, we 
obtain that 2oq ^ at, which concludes in the second case. 

Third case, we assume only that 7 > f (and always /3 ^ f). Let h £ §(77) be the 
unique point on the spherical segment If, f 1 such that d§(q,h ) = 7 In particular, we have 
that *(/, f) = ds(/, h) + d§{h, /') and' 

at = (k{ft, ft) < <k{ft, h t ) + d§(h t ,f), 

def 

where h t = T(g, h,t). The first case shows that d§(f t ,h t ) ^ d§(f,h) and the second case 
shows that d§(h t , f' t ) ^ 2 d§(h, f). Thus a t < d§(/, h) + 2 d§(h, f) < 2 d§(f, f). □ 


Recall that M(f,g,Q denotes sup teJ g t , see §1.2 and Equation (4). 

Lemma 16. If D 3 / 2 M(f, g, Cf) 2 p ^ then J = J' = [0,1] and for any t £ [0,1]: 

(i) (1 + £)~V t < IH < (1 + OmD 

(ii) D 3 / 2 g, t dp(f t , CD < gy- 


Proof. The assumption implies that M{f,g,Cf) < 00 , and thus J = [0,1], by Lemma 4. Let S' 
the set of all t £ J' such that D 3 ^ 2 p t dp(Ct, Ct) ^ si- H is a nonempty closed subset of J'. 
Let t £ S. By Lemma 15, we have dp(f t , f' t ) ^ 2 p, so 


D 3 ^ 2 Ptdp(ft, ft) < 


£ 

4(1 + e)' 


Proposition 5 implies that there exists a root 77 of f' t such that dp(g, ( t ) ^ 2(1 + e)pt.p and 
(1 + £)“ V* < M/D v) < (1 + e)lH- Because d P (?7, CD < d P (?7, Ct) + M»(Ct, CD and t £ S we 
obtain 


D 3 ^ 2 p(f,v)dp(g, CD < ZM 2 (l+£)/i t ^2(1 + £)/x t p + ^7)372)7^ ^ (I+ e ) 2 y^+(1+ £ ) <1 -• 


and Theorem 2 implies that Ct approximates 77 as a root of fj.. Since it is also an exact 
root of /D this implies Ct = V- In particular D 3 ^ 2 p t dp{Ct,Ct) ^ 2(1 + e)D 3 l 2 pf p < H ^ . 
Thus f is in the interior of S, which proves that S is open and finally that S = J. Moreover, 
since M ^ (1 + e)pt, M I s bounded on J’, thus J’ = [0,1]. □ 


This leads to the procedure HC / , see Algorithm 2. It modifies procedure HC (Algorithm 1) 
in only one respect: each iteration checks up on the failure condition D 3 / 2 p(h, z ) 2 p>y. If 
the failure condition is never met, then HC / computes exactly the same thing as HC. Recall 
that M(f, g,rf) denotes the maximum condition number p that arises in the homotopy 
continuation HC(//< 7 , 77 ), and that I p (f,g, rf) denote the integral of pP along the homotopy 
path from g to /, see §1.2 and Equation (4). 
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Algorithm 2. Homotopy continuation with precision check 
Input. f,g£ §("%), zeP“ and p > 0. 

Output, tu £ P" or FAIL. 

Specifications. See Proposition 17. 

function HC'(/, g, z, p) 

t £- 1/ (101 L> 3/2 /i(g, z) 2 d§(f, g)) 
h <— g 

while 1 > t and D 3 I 2 p(h, z) 2 p ^ Ay do 
h T(g, f,t) 

z AT(h, z ) 

t <r- t + 1/ (l01L> 3/2 /x(/i, z) 2 d$(f, g)) 

end while 

if D 3 / 2 p(h, z) 2 p > Ay then return FAIL 
else return 2 
end if 

end function 


Proposition 17. If ds(f, g) < \ and d(f,f') ^ p, then the procedure HC '(f',g,g,p): 

(i) terminates and performs at most 158 D 3 I 2 d§(f : g)/ 2 (/, g, v) +4 steps; 

(ii) outputs an approximate root of f, or fails; 

(Hi) succeeds (i.e. outputs some z £ P™ ) if and only if D 3 I 2 M(f ,g,g) 2 p ^ Ay,- 
(iv) succeeds if D 3 I 2 M(f 1 g,g) 2 p ^ Ag. 

Proof. At each iteration, the value of t increases by at least 151p/ (101cZs( < /’ / , <7)), thus there 
are at most 101rfs(y / , <7)/(151/a) iterations before termination. 

By construction, the procedure KC'(f', g, 77, p) fails if and only if at some point of the 
procedure HC(/', g, 77, p) it happens that D 3 ! 2 p(h, z) 2 p > y|y, In other words, the proce¬ 
dure HC 1 (f, g,g, p) fails if and only if D 3 ^ 2 M{f , g,g) 2 p > yyy, by definition of M. And 
since the procedure terminates, it succeeds if and only if it does not fail. This proves (iii). 

Let us bound the number K'(f , g , g, p) of iterations of the procedure HC '(f, g , g, p) before 
termination. If HC^/', g, g, p) succeeds, then K'(f, <7, g, p) = K(f,g,g). Furthermore 

K'{f',g,g,p) = sup{K(f' s ,g,g) \ s£ [0,1], HC'(/(, g, g, p) succeeds} . (16) 

Let s € [0,1] such that HC '(f' s ,g,g,p) succeeds, that is to say D 3 I 2 M(f' s , g, g) 2 p < Ay. 
Theorem 8(ii) shows that 

(1 + £)~ 2 M(f' s ,g,g) ^ M(f s ,g,g) < (1 + e) 2 M{f' s , g, g). 

In particular D 3 I 2 M(f s ,g,g) 2 p ^ A±£l_ ^ ypj and Lemma 16 shows that (1 + e)~ 2 < p! t ^ 
(1 +e)p t for all t ^ s. So we obtain that (l+e)“ 2 / 2 (/ s , g, g) ^ h{f s ,g,ri) < (l+e) 2 / 2 (/ s , g, g) 
and 

K(fs,9,v) < 136D 3/2 d§(f' s ,g)I 2 (f' s ,g,ri) by Theorem 8(iii) 
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136(1 + e) 2 D 3/2 (d§(f s ,g) + 2p) I 2 (fs,g,g) by Lemma 15. 

Besides D 3 / 2 I 2 (f s ,g,g)p ^ (1 + e) 2 D 3 ^ 2 M(f' s , g,g) 2 p ^ B+D , so we obtain 

K(fs,g,v) l58D 3/2 d§{f s ,g)I 2 (f s ,g 1 g) +4 < 158D 3/2 d§(f 1 g)I 2 (f,g,g) +4. 

Together with Equation (16), this completes the proof of (i). 

Let us assume that the procedure g, g, p) succeeds and let 2 be its output, which 

is nothing but HC(/', g, 77 ). Theorem 8 (v) shows that D 3 ^ 2 fi[dp(z, £() < where is 
the root of /( = /' obtained by homotopy continuation. As above, with s = 1, we check 
that /xi < (1 + c)/i'| and H 3 / 2 /x' 1 dp(Ci,Ci) ^ gy using Lemma 16. Thus 

D 3 ^ 2 gidp(z, Ci) ^ (1 + e) ^ < B -. 

Then z approximates Ci as a root of / 1 , by Theorem 2. This proves (ii). 

Lastly, let us assume that D 3 / 2 M(f,g,rj) 2 p ^ Lemma 16 implies that M(f,g,rj) ^ 

(1 + g, 77 ) and Theorem 8 (ii) shows that M(f,g,i 7 ) < (1 + e) 2 M(f , g,r]). Thus 

D 3 / 2 M(f,g, V ) 2 p^ (l+s) 6 D 3 / 2 M(f,g,g) 2 p^ <Bj J_ 

and HC'(/', g, 77 , 77 ) succeeds. This proves (iv). □ 

2.3 A deterministic algorithm 

Let / £ S("H) be the input system to be solved and let Q ^ 1 be a given precision. We 
compute 

/' = L/J Q, (g,v) = BP({/} Q ), e = sign( 7 r /2 - d§(f,g)) and p = 3N 1/2 /Q. 

Lemma 11 shows that dg(/, f) ^ p. Then we run the homotopy continuation procedure 
with precision check eg, 77 , p), which may fail or output a point z £ P n . If it does 

succeed, then Proposition 17 ensures that z is an approximate root of /. If the homotopy 
continuation fails, then we replace Q by Q 2 and we start again, until the call to HC 7 succeeds. 

This leads to the deterministic procedure DBP, Algorithm 3. If the computation of DBP(/) 
terminates then the result is an approximate root of /. Section 2.4 studies the average 
number of homotopy steps performed by DBP(/) while Section 2.5 studies the average total 
cost of an implementation of DBP in the BSS model extended with the square root. 

2.4 Average analysis 

Let f £ S (T~L) be the input system, a uniform random variable, and we consider a run of the 
procedure DBP(/). Let Q k be the precision at the fcth iteration, namely Q k = N 2 . We set Q k 

also fk, gk, Vk, £fc, 

fk = L f\Qk . ( gk,Vk ) = BP({/}qJ, £fc = sign( 7 r /2 - d s (f,g k )) and p k = 3 N 1/2 /Q k . Pk 

Let Q be the least k such that the homotopy continuation with precision check HC ' {fk, £kgk, 77 k, Pkfl 
succeeds. Note that Q is a random variable. To perform the average analysis of the total 
number of homotopy steps, we first deal with each iteration separately (Lemmas 18 and 19) 
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Algorithm 3. Deterministic variant of Beltran-Pardo algorithm 

Input. / G H 
Output. z6P" 

Postcondition, z is an approximate root of / 

function DBP(/) 

Q^N 

repeat 

Q t— Q 2 
f L/J Q 

£ <- sign(Re(/, g)) 

P 4r- ( 2 Nf^/Q 

z <- HC '(f,eg,g,p) 
until HC / succeeds 

return 2 
end function 


and then give tail bounds on the probability distribution of Q (Proposition 20). Even if the 
number of steps in each iteration are not independent from each other and from 0, Holder’s 
inequality allows obtaining a bound on the total number of steps (Theorem 21). 

Let (g,rf) G V he a random variable with distribution p st d and independent of /. g, g 

Lemma 18. Let 0 : % x V —> K. be any nonnegative measurable function. For any k ^ 1, 

E (0(/fc,£fc3fe,r?fc)) < 10E (0(/fc, 5 , 77 )). 

Proof. It is an application of Lemma 11. We first remark that £ k G {—1,1} so 
6 ( fk,£k9k, Vk) < © (fk, 9k, Vk) + 0 (fk, -9k, 9k) ■ 

Then 

E(e(} M )) = j sm) e (l/j 0 ., bp ({ /[,. )) if 

exp ( 2 q V ~ ) r 

< — fkjfTyd / 0(L/jQ fc > BP (5))d/d3 by Lemma 11 

vol(S (H)) 2 J S («)xS(«) 

= vol(§{H))^ / Jv 6 *- L^J Qk ’ 9 ' ^ d ^ df>std ^ 9, ^ by Theorem 9 
= ex P E ( 0 (/fc’ 9, v)) ■ 

Similarly, E (©(/*,, ~9k, Vk)) exp E (0(/ fc , -g, rff), and since g and -g have the 

same probability distribution, E (0(/fc, — g, g)) = E (0(/fc, g, g)). To conclude, we remark 
that Q k ^ TV 2 and that ^5. □ 

Lemma 19. K(I p (f,g,g)) = E (g(g,g) p ) for any p ^ 1 and k^l. 
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Proof. Let ht. = T(g,f,t), for t £ [0,1], and let ft be the associated homotopy contin¬ 
uation. Let t £ [0,1] be a uniform random variable independent from / and ( g,g ). 

Clearly E(I p (f, g, g)) = E (p(h T , f T ) p ), so it is enough to prove that (h T ,fr) ~ Pstd- The 
systems f and g are independent and uniformly distributed on 8>(fH). So their probability 
distributions is invariant under any unitary transformation of PL. Then so is the probability 
distribution of h t for any t £ [ 0 , 1 ], and there is a unique such probability distribution: the 
uniform distribution on S (PL). The homotopy continuation makes a bisection between the 
roots of g and those of h t . Since g is uniformly chosen among the roots of g, so is ft among 
the roots of h t . That is, ( h t ,ft ) ~ Pstd for all t £ [0,1], and then ( h T ,f T ) ~ p s td- □ 

Proposition 20. P(fi > k) ^ 2 1 ' D 9 ^ 4 n 3 ^ 2 N 7 ^ 4 Q k 1 ^ 2 . 

Proof. The probability that fl > k is no more than the probability that HC / (f k , g k ,gk, Pk) 
fails. By Lemma 18, P (KG'{f k ,e k g k ,g k , p k ) fails) < 10P (UC'(f k ,g,g,p k ) fails). Given 
that ds(f,f k ) < Pk, 

P (HC '(fk,9,V,Pk) fails) < P ^D 3/2 M(f,g,g) 2 p k Js by Proposition 17(iv) 

< P ^ D 9/2 I 3 (f , g, g) 2 p k ^ 236 7 ^ 2 ) by Proposition 7 

^ 151 • 236 1 / 2 D 9 / 4 p 1 J 2 E (/ 3 (/, g , 77 )) by Markov’s inequality. 

Lemma 19 and Theorem 10 imply then 

E [W, 9, v)) < E (p(g, g) 3 ) < 3(nlV) 3 / 2 . 

All in all, and since p k = 3 N 1 / 2 /Q k , 

P(fl > k) < 10-(151-236 1/2 D 9/4 )-(31V 1/2 /Q fc ) 1/2 -3(nA') 3/2 2 47 D 9/4 n 3/2 N 7/4 Qf 1/2 □ 

Let K(f) be the total number of homotopy steps performed by procedure DBP(/) and let K{f) 
the number of homotopy steps performed by procedure HC ' (fk,£k9k,Vk, Pk) be denoted by 
K'(fk,£kgk,Vk,Pk), so that 

n 

K(f) = £ K'{fk,£kgk,Vk, Pk), 
k= 1 

Theorem 21. If N ^ 21 then E(K(f)) £, 2 17 nD 3 / 2 N. 

Proof. Let X k = K'(f k ,s k g k ,g k , p k ) and let 0 < p < |. By Lemma 18 and Proposition 17(i), 

« \ p\ 1 /p 

158 D 3 / 2 d$(f,g)I 2 (f,g,r 1 ) +4) J , 
and because d§(f,g) ^ 7 r and by Minkowski’s inequality, we obtain 

< 10 (l58L> 3/2 7rE(/ 2 (/, 5 ,77) p ) 1/p + 4) . 

Jensen’s inequality implies that I 2 (/, g, g) p < I 2p {f, g, g). Then E {I 2p {f, g, g)) < A= 2 p{ nN Y < 
3(nN) p , by Lemma 19 and Theorem 10. In the end, 

E(A£) 1/p 15000 nD 3/2 N. (17) 


20 



Now, let p = log TV/(log N — 1). If N ^ 21 then p |. We write the expectation of K (/) as 

17 oo 

=E(J2Xk) = E 

h—1 k -1 

Let q = 1/log-ZV, so that ^ + | = 1. From Holder’s inequality, E(Afcln^fc) ^ E(X^) 1 / p P(fl ^ 
k) x l q and thus 

OO 

E(tf(/)) < maxE (X^E P(fi ^ fc) 1/9 . 

^ fc=l 

Lemma 22 below, with C = 1, L = 4 and d = l/q, shows that 
00 _ 9 17 / lo s 21 p 5 

^ fc) 1 / 9 ^ L + 1 H-^6- 

ti e ~ 1 

The claim follows then from Equation (17) and 6 • 15000 2 17 . □ 


Lemma 22. For any C, 5 > 0 and any integer L ^ 2 such that C < N s2 , 

oo l+ i (^n r ,L + 2 

C k P(n > k) s < C k + ' 

fe=1 A—1 


JV 521, - c 


Proof. For any k, Proposition 20 implies that 

P(fi > k) = P(fi > k - 1) < 2 17 D 9/4 n 3/2 AT 7/4 Q“^{ 2 < 2 17 N 5 N~ 2k ~\ 

using D ^ N and ?r 2 ^ AT. Moreover, 2 P_1 ^ p, for any integer p, so that N~ 2k 2 ^ 
N~ 2 (fc-^- 1 ). Of course, it also holds that P(f2 ^ k) ^ 1. Thus 

oo L+l oo 

C fc P(H > fc) 5 < C ' k + (2 17 iV 5 )' 5 C k N -82 L (k-L-l), 

fc^l fc=l k=L+2 

and the latter sum is a geometric sum which evaluates to C L+2 / (N S2h — C). □ 


2.5 Implementation in the BSS model with square root 

Algorithms HC / and DBP (Algorithms 2 and 3 respectively) have been described assuming 
the possibility to compute exactly certain nonrational functions: the square root, the 
trigonometric functions sine and cosine and the operator norm of a linear map. A BSS machine 
can only approximate them, but it can do it efficiently. I propose here an implementation 
in the BSS model extended with the ability of computing the square root of a positive real 
number at unit cost. We could reduce further to the plain BSS model at the cost of some 
lengthy and nearly irrelevant technical argumentation. We now prove the main result of this 
article: 

Theorem 23. There exists a BSS machine A with square root and a constant c > 0 such 
that for any positive integer n and any positive integers d\, ..., d n : 

(i) A(f) computes an approximate root of f for almost all f £ TL; 
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(ii) if f £ S(H) is a uniform random variable, then the average number of operations 
performed by A(f) is at most cnD 3 / 2 N(N + n 3 ). 

Firstly, we describe an implementation of Algorithms HC / and DBP in the extended BSS 
model. The first difficulty is the condition number p(f,z): it rests upon the operator norm 
for the Euclidean distance which is not computable with rational operations. While there 
are efficient numerical algorithms to compute such an operator norm in practice, it is not 
so easy to give an algorithm that approximates it in good complexity in the BSS model. 26 
Fortunately, we can easily compute the operator norm of a matrix M £ C" x ™ within a 
factor 2 as follows: 2 ' we first compute a tridiagonalization T of the Hermitian matrix M t M 
with 0(n 3 ) operations, using Householder’s reduction, and then 28 

where ||T||i is the operator fi-nornr of T, that is the maximum fi-norm of a column. Therefore, 
up to a few modifications in the constants, we may assume that p(f, z ) is computable in 0(n 3 ) 
operations, given d f(z). 

The second difficulty lies in the use of the trigonometric functions sine and cosine. They first 
appear in the definition of the geodesic path T, Equation (3), which is used in Algorithm 2. 
In the case where d§(f,g) < tt/2, it is good enough to replace T(g, /, S) by 

5f + (l-6)g 

n*/+(i-%ir 

This is classical and implies modifications in the constants only. 29 The trigonometric functions 
also appear in Sibuya’s function S, see Equation (11). This issue can be handled with power 
series approximations: 

Lemma 24. There is a BSS machine with square root that computes, for any N and 
any x £ [0, l] 2Ar_1 , with 0(N logN) operations, a point S(x) £ E> 2N ~ 1 such that 

[ Q{S(x))dx < u ^ 2 N-i\ [ 0 (y) d 2b 

J[ 0,1]2« —1 VOl(S 2JV J§2N-l 

for any nonnegative measurable function 0 : S 2W_1 —> R. 

Sketch of the proof. For any positive integer Q , let Fq(x) be the Taylor series expansion 
at 0, truncated at x®, of the entire function (exp(2z7rx) — l)/(x — 1). It is a polynomial 
of degree Q that can be computed with O(Q) operations, assuming that tt is a constant 
of the machine, by using the linear recurrence (n + 2 )u n +2 = (2z7t + n + 2)u„+i + 2 ; mu n 
satisfied by the coefficients of Fq. Let Cosg(a;) and Sing (a:) be the real and imaginary parts 
of (1 + (x - 1)Fq(x))/|(1 + (x - l)Fg(a;))| respectively. 

The function x £ [0,1] —> (Cosq(x), Sing (a;)) gives a parametrization of the circle S 1 whose 


26 See for example Armentano, Beltran, Biirgisser, Cucker, and Shub, A stable, polynomial-time algorithm for 
the eigenpair problem or Armentano and Cucker, “A randomized homotopy for the Hermitian eigenpair 
problem”; unfortunately the Gaussian distribution that they assume does not fit the situation here. 

27 I thank one of the referees for having communicated this method to me. 

28 Kahan, Accurate eigenvalues of a symmetric tri-diagonal matrix. 

29 See for example Biirgisser and Cucker, Condition , §17.1. 
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Jacobian is almost constant: we can check that there is a universal constant C > 0 such that 


|Cosg(x) 2 + Sing (a :) 2 — (27r) 2 | ^ Ce ®. 

Thus for any continuous function 6 : S 1 —> R. 

r 1 1 + Ce~ Q f 

/ 0(Cosg(x),Sing(x))dx < --- / 0(y)dy. 

Jo 2 tt J § i 

Let S be the function [0, 1] 2JV_1 —>. S 2W_1 defined in the same way as S, Equation (11), 
but with Cosg and Sing in place of sin and cos respectively, with some Q ~ log IV such 
that (1 + Ce~®) N ^ 2. It is easy to check that S satisfies the desired properties. □ 

In Algorithm DBP, there is no harm in using S in place of S. We obtain this way variants 
of Algorithms HC / and DBP that fit in the BSS model with square root. It remains to 
evaluate the overall number of operations. It is well known that f(z) and d f(z) can be 
computed at a point z £ C " +1 in O(N) operations—the latter as a consequence of a theorem 
of Baur and Strassen. Together with the approximate computation of the operator norm 
discussed above, this implies the following: 

Lemma 25. There exists a BSS machine with square root that compute n(f,z) (within a 
factor 2) and J\f(f , z), for any f £ TL and z £ P", in 0(N + n 3 ) operations. 

The cost of the fcth iteration in Algorithm DBP is dominated by the cost of computing |_/J Q k 
and BP({/}g fc ) and the cost of the call to HC’. The cost of the call to HC’ is dominated 
by the cost of the homotopy steps. Each homotopy step costs 0(N + n 3 ) operations, by 
Lemma 25. 

We now evaluate the cost of computing [/J<2 an d BP({/}g). Naturally, the integral 
part A(x) of a real number x is not a rational function of x but it can be computed in the 
BSS model in 0(log(l + |x|)) operations using the recursive formula, say for x ^ 0, 


TO if x < 1 

A(x) = l 2A(x/2) if x < 2A{x/2) + 1 
l_2A(x/2) + l otherwise. 

Hence we can compute Q~ x A(Qx) in O(logQ) operations, for any positive integer Q and x £ 
[0,1]. It follows that one can compute [/Jq i 11 O(NlogQ) operations. The computation 
of {/}g is similar and it is done in 0(N log Q + NlogN) operations, by Lemma 24, using S 
in place of Sibuya’s function S. Finally, given {/}g, one compute BP({/}g) in 0{N 2 ) 
operations. In the end, the cost of the /cth operation, excluding the call to HC 7 , is thus 0(N 2 + 
N log Qk) operations. 

Hence, the overall cost of the algorithm is 

O ^(N + n 3 )K(f) + J2( N2 + N ^°sQk) ] j =o(^N + n 3 )K(f) + N 2 n + NY^ogQ^j , 

where K(f) is the total number of homotopy steps. By Theorem 21, E(A'(/)) = 0(nD 3 / 2 N), 
so it only remains to bound the expectations of fi and l°gQfc- 

30 Baur and Strassen, “The complexity of partial derivatives”. 
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Lemma 26. E(fi) ^ 7 

Proof. By Lemma 22, with C = 1, 6 = 1 and L = 5, 

o L+17 ivr5 

^ ^ k) < L + 1 + 2I , _ 7' □ 

fc=i 

Lemma 27. E ^X^fc=i l°g Qk) ^ 129 log TV. 

Proof. Because Qk = N 2k , 

Q oo oo 

E( log Qk) = J2 lQ g Q k p (^ > k ) = log TV ^ 2 fc P(fi > k). 

k=i fe =i fc=i 

Lemma 22, with <5 = 1, C = 2 and L = 5 gives that 

00 oL+19/\r5 

5Z 2fc ^ k ) < 2 i+2 + ^ _ 2 129, 

fc=i 

where we used that TV ^ 2. □ 

This concludes the proof of Theorem 23. 

Conclusion 

The derandomization proposed here relies on extracting randomness from the input itself, 
which is made possible by the BSS model and the infinite precision it provides. Actually, this 
might also work under finite, and rather moderate, precision. Indeed, in the fcth iteration 
of Algorithm 3, we need, very loosely speaking, about logQfc digits of precision but not 
much more, because by construction, the homotopy continuation procedure HC / aborts when 
more precision would be required for the result to be relevant. And then, Lemma 27 shows 
that logQfc, is typically no more that 129 log TV. So it is reasonnable to think that, extending 
the work of Briquel et al., 31 we may run a variant of Algorithm 3 on a finite precision machine 
and obtain a significant probability of success as long as we work with CTogTV digits of 
precision, for some constant C > 0. 

Besides, Armentano et al. 32 proposed recently a new complexity analysis of Beltran-Pardo 
algorithm which relies on a refined homotopy continuation algorithm. 33 They obtained a 
randomized algorithm that terminates on the average on a random input after 0(nD 3 ^ 2 N 1 ^ 2 ) 
homotopy steps. This is a significant improvement on the previously known 0(nD 3 ^ 2 N) 
bound. The derandomization method should also apply to this refined algorithm, in all 
likelihood, but this is not immediate: to devise the homotopy continuation with precision 
check, we had to look deep inside the continuation process. Adapting the method to a refined 
homotopy continuation process will inevitably lead to further difficulties. 

31 Briquel, Cucker, Pena, and Roshchina, “Fast computation of zeros of polynomial systems with bounded 
degree under finite-precision”. 

32 Armentano, Beltran, Biirgisser, Cucker, and Shub, “Condition length and complexity for the solution of 
polynomial systems”. 

33 Shub, “Complexity of Bezout’s theorem. VI. Geodesics in the condition (number) metric”; Beltran, “A 
continuation method to solve polynomial systems and its complexity”. 
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